% rc9 Feb-2011 - 
%  data quality checks - b_norm, e_zero, esum, e_avg & e_03d
%
% writes a ".gcq" (good channel quality) file (.gcq + quality selection)
%  for use by mtg_good.m

function e_sum_qual(in_path,out_path,log_file)

global hdr file_name i_block dc ac ph adc aac aph n2 ssdc ssac ssph
global nir_ph nir_ac nir_dc msg zap_flg id_tags evt_PathName avg_PathName
% global good_dc_pts good_dc_chs good_ac_pts good_ac_chs
% global good_ph_pts good_ph_chs event n_good
% good_chs

clear out_path

fprintf('e_sum_qual\n');

print_flg=0; % 0=> 'don't print'; 1=> 'print to screen'; 2=> 'print to file';
if print_flg==1
    pfid=1;
elseif print_flg==2
    pfid=fopen([avg_PathName 'qual_log_' file_name],'w');
end

% in_path=[hdr.data_path hdr.exp_name hdr.data_type hdr.hi '-' hdr.lo '\'];
% out_path=[hdr.data_path hdr.exp_name hdr.data_type hdr.hi '-' hdr.lo '\'];

e_sum_id = '$Id: e_sum.m 488 2005-12-15 20:15:45Z edmaclin $';
check = strcmp(e_sum_id,id_tags); % Make sure we don't already have the tag
if sum(check) == 0
    id_tags{size(id_tags,1)+1,1} = e_sum_id;  % Add it in then.
end

% hdr.n_good.tot=zeros(1,hdr.n_chans);

scrsz = get(0,'ScreenSize'); % for debug plotting

% Determine whether input and output files exist
%  !!! 'exist' function doesn't work in compiler :-(
%
disp(['in_path = ' in_path file_name])
in_file_flg=fopen([in_path file_name]);
evt_file_name=[evt_PathName '\' file_name];
disp(['Event file = ' evt_file_name])
evt_file_flg=fopen([evt_file_name]);
k=findstr('.',file_name);
out_file_name=[file_name(1:k(1)) 'avm'];
nir_file_name=[file_name(1:k(1)) 'nir'];
xls_file_name=[file_name(1:k(1)) 'xlsx'];
gcq_file_name=[file_name(1:k(1)) 'gcq'];

out_file_flg=fopen([avg_PathName out_file_name]);

if in_file_flg > 1, fclose(in_file_flg); end
if evt_file_flg > 1, fclose(evt_file_flg); end
if out_file_flg > 1, fclose(out_file_flg); end

if ~strcmp(log_file,'')
    logid = fopen(log_file,'a');
else logid = 1; end

fprintf(logid,'Summing...\n');

if out_file_flg < 0
	if ((in_file_flg > 0) & (evt_file_flg > 0))       
        % Load 'event' matrix in .mat format 
        %
        msg=sprintf('Reading %s',evt_file_name);
        load(evt_file_name,'event','-mat');
        fprintf(logid,'Reading %s\n',evt_file_name);

        % load EROS data in .mat format
        %
        temp=[in_path file_name];
        load(temp,'-mat');
        fprintf(logid,'Reading %s\n',temp);

        % Sum entire block (for hemodynamic & spectral analyses)
        % NOTE: May remove or make optional.  Not usually used.
%         fprintf(logid,'Summing blocks...\n');
        nir_ph=nir_ph+ph;
        nir_ac=nir_ac+ac;
        nir_dc=nir_dc+dc;
        
        n_trials=size(event,1);
        n_chans=size(dc,2);
        
        % quality check -
        %
        % check each channel for BLOCK mrac, mrdc
        % done here (in addition to b_norm) to allow use of liberal "norming"
        %
        % disp('setting block good_mrdc&ac chs')
        
        % first reapply norm thresholds (in case avg threshs are liberal (0!)
        %
        good_norm_mrdc_chs=find(mrdc>hdr.qual.norm_mrdc_thr);
        good_norm_mrdc_flgs=zeros(1,hdr.n_chans); % make list for .xls file
        good_norm_mrdc_flgs(good_norm_mrdc_chs)=1;
        
        good_norm_mrac_chs=find(mrac>hdr.qual.norm_mrac_thr);
        good_norm_mrac_flgs=zeros(1,hdr.n_chans); % make list for .xls file
        good_norm_mrac_flgs(good_norm_mrac_chs)=1;
        
        % now apply avg dc & ac thresholds
        %
        good_mrac_chs=find(mrac>hdr.qual.mrac_thr);
        
        hdr.n_good.mrac(good_mrac_chs)=hdr.n_good.mrac(good_mrac_chs)+1;
%         fprintf('%3.0f good_mrac_chs\n',length(good_mrac_chs));
%         fprintf('%3.0f good_chs\n',length(good_chs));
        good_mrac_flgs=zeros(1,hdr.n_chans); % make list for .xls file
        good_mrac_flgs(good_mrac_chs)=1;
        
        good_mrdc_chs=find(mrdc>hdr.qual.mrdc_thr);
        
        hdr.n_good.mrdc(good_mrdc_chs)=hdr.n_good.mrdc(good_mrdc_chs)+1;
%         fprintf('%3.0f good_mrdc_chs\n',length(good_mrdc_chs));
%         fprintf('%3.0f good_chs\n',length(good_chs));
        good_mrdc_flgs=zeros(1,hdr.n_chans); % make list for .xls file
        good_mrdc_flgs(good_mrdc_chs)=1;
        
%         if print_flg>0
%         fprintf(pfid,'%3.0f / %3.0f good_mrdc_chs\n',length(good_mrdc_chs),n_chans);
%         fprintf(pfid,'%3.0f / %3.0f good_mrac_chs\n',length(good_mrac_chs),n_chans);
%         end

        % Compute sums per condition and channel
        %
        % disp('looping through channels')
        for i_trial=1:n_trials
            
            % re-initialize "good_chs" for each trial
            %
            % First select by distance (all chs if not checked in gui)
            % "good_dist_chs" must exist in "norm" file ...
            %
            good_chs=good_dist_chs;
%             fprintf('%3.0f good_dist_chs\n',length(good_dist_chs));
    
            good_chs=intersect(good_chs,good_norm_mrdc_chs);
            good_chs=intersect(good_chs,good_norm_mrac_chs);
            good_chs=intersect(good_chs,good_mrac_chs);
            good_chs=intersect(good_chs,good_mrdc_chs);
%             fprintf('%3.0f good_chs\n',length(good_chs));

            i_cond=event(i_trial,2); % only include trials for current condition
            start=event(i_trial,1)-hdr.avg_baseline+1;
            stop=event(i_trial,1)+hdr.epoch_length;

            if (start<1) | (stop>length(dc))
%                 fprintf('e_sum_qual 198; i_trial=%4.0f,  start=%4.0f,  stop=%4.0f, length(dc)=%4.0f\n',...
%                     i_trial,start,stop,length(dc))
            else
                % quality control:
                %
                % check each channel for dc range
                %
                % disp('Setting good trial chans')
                dc_max=max(dc(start:stop,:));
                dc_min=min(dc(start:stop,:));
                dc_range=dc_max-dc_min;
                good_dc_range_chs=find(dc_range<hdr.qual.dc_range);
                good_chs=intersect(good_chs,good_dc_range_chs);
                hdr.n_good.dc_range(good_dc_range_chs)=hdr.n_good.dc_range(good_dc_range_chs)+1;
%                 fprintf('%3.0f good_dc_range_chs\n',length(good_dc_range_chs));
%                 fprintf('%3.0f good_chs\n',length(good_chs));
                
                % check each channel for ac min
                %
                % disp('checking AC min')
                ac_epc_min=min(ac(start:stop,:)+1); % restore mean = 1
                good_ac_chs=find(ac_epc_min>=hdr.qual.ac_min);
                good_chs=intersect(good_chs,good_ac_chs);
                hdr.n_good.ac(good_ac_chs)=hdr.n_good.ac(good_ac_chs)+1;
%                 fprintf('%3.0f good_ac_chs\n',length(good_ac_chs));
%                 fprintf('%3.0f good_chs\n',length(good_chs));
                
                % check each channel for delay range
                %
                % disp('checking delay range')
                ph_max=max(ph(start:stop,:));
                ph_min=min(ph(start:stop,:));
                ph_range=ph_max-ph_min;
                good_ph_chs=find(ph_range<hdr.qual.delay_range);
                good_chs=intersect(good_chs,good_ph_chs);
                hdr.n_good.ph(good_ph_chs)=hdr.n_good.ph(good_ph_chs)+1;
%                 fprintf('%3.0f good_ph_chs\n',length(good_ph_chs));
%                 fprintf('%3.0f good_chs\n',length(good_chs));
                
                % check each channel for modulation min for this epoch
                % low modulation may mean AC clipping
                % 1st compute modulation for this epoch
                % mod is meaningless (-, NaN) if there is little light
                % so only check channels with dc>hdr.qual.mrdc_thr, OR 
                % > 100 if hdr.qual.mrdc_thr<100 (mainly for testing)
                %
                tmrdc=repmat(mrdc,stop-start+1,1); % mrdc for each pt in epoch
                
                % fix 0 or negative DC
                %
                if hdr.qual.mrdc_thr<100
                    tmp=find(tmrdc<100); 
                    tmrdc(tmp)=100;
                else
                    tmp=find(tmrdc<hdr.qual.mrdc_thr);
                    tmrdc(tmp)=hdr.qual.mrdc_thr; 
                end
                
                tmrac=repmat(mrac,stop-start+1,1); 
                tmrac(tmp)=0;
                epc_mod=(tmrac.*(ac(start:stop,:)+1))./(tmrdc.*(dc(start:stop,:)+1));
                %
                % use 'min' or 'mean' modulation during epoch ?
                %
%                 min_mod=min(epc_mod); % minimum mod on each channel
%                 good_mod_chs=find(min_mod>hdr.qual.mod_thr);
                mean_mod=mean(epc_mod); % mean mod on each channel
                tmp=find(isnan(mean_mod));
                mean_mod(tmp)=0; % clean up remaining NaNs... ???
                good_mod_chs=find(mean_mod>=hdr.qual.mod_thr);
                
%                 if print_flg>0
%                     fprintf(pfid,'\nmean_mod\n');
%                     fprintf(pfid,'%5.2f ',mean_mod);
%                     fprintf(pfid,'\ngood_mod_chs\n');
%                     fprintf(pfid,'%3.0f ',good_mod_chs);
%                     fprintf(pfid,'\n');
%                     neg_mod_chs=find(mean_mod<0);
%                     mean_mod(neg_mod_chs)=0;
%                     bad_mod_chs=find(mean_mod<hdr.qual.mod_thr);
%                     fprintf(pfid,'bad_mod_chs = %3.0f, ',bad_mod_chs);
%                     fprintf(pfid,'\nmean_mod(bad_mod_chs:\n');
%                     for i_tmp=1:length(bad_mod_chs)
%                         fprintf(pfid,'%3.0f) %5.0f\n',bad_mod_chs(i_tmp),mean_mod(bad_mod_chs(i_tmp)));
%                     end
%                 end
                
                good_chs=intersect(good_chs,good_mod_chs);
                hdr.n_good.mod(good_mod_chs)=hdr.n_good.mod(good_mod_chs)+1;
%                 fprintf('%3.0f good_mod_chs\n',length(good_mod_chs));
%                 fprintf('%3.0f good_chs\n',length(good_chs));
                
                % increment accepted channels for this bin
                %
                n2(good_chs,i_cond)=n2(good_chs,i_cond)+1;
                %
                % and total
                %
                hdr.n_good.tot(good_chs)=hdr.n_good.tot(good_chs)+1;

                % sum good channels for this trial
                %
                adc(:,good_chs,i_cond)=adc(:,good_chs,i_cond)+dc(start:stop,good_chs);
                ssdc(:,good_chs,i_cond)=ssdc(:,good_chs,i_cond)+dc(start:stop,good_chs).^2;
                aac(:,good_chs,i_cond)=aac(:,good_chs,i_cond)+ac(start:stop,good_chs);
                ssac(:,good_chs,i_cond)=ssac(:,good_chs,i_cond)+ac(start:stop,good_chs).^2;          
                aph(:,good_chs,i_cond)=aph(:,good_chs,i_cond)+ph(start:stop,good_chs);
                ssph(:,good_chs,i_cond)=ssph(:,good_chs,i_cond)+ph(start:stop,good_chs).^2;   

            end % epoch within data block
            
            % print running total of accepted trials/ch
            %
%             current_good_chs=zeros(1,hdr.n_chans);
%             current_good_chs(good_chs)=1;
%             fprintf('%3.0f)',i_trial); 
%             fprintf('%3.0f ',current_good_chs);
%             fprintf('\n    ');
%             fprintf('%3.0f ',hdr.n_good.tot);
%             fprintf('\n');

        end % i_trial
       
        % create .gcq file reflecting dist & quality selections
        %
        % N.B. this is a list of the number of accepted epochs in all
        % channels, NOT a list of "good" channels !!!
        %
        n_good_tot=hdr.n_good.tot;
%         keyboard
        disp(['Saving ' hdr.data_path hdr.exp_name 'reg\' gcq_file_name])
%         save([hdr.data_path hdr.exp_name 'reg\' gcq_file_name],'n_good_tot','-ascii')
        fid=fopen([hdr.data_path hdr.exp_name 'reg\' gcq_file_name],'w');
        fprintf(fid,'%4.0f\n',n_good_tot);
        fclose(fid);
        
        if print_flg>0
            fprintf(pfid,'Norm parmaters:\n');
            fprintf(pfid,'hdr.qual.ph_out_thr= %5.2f\n',hdr.qual.ph_out_thr);
            fprintf(pfid,'hdr.qual.norm_mrdc_thr= %5.2f\n',hdr.qual.dc_max);
            fprintf(pfid,'hdr.qual.norm_mrac_thr= %5.2f\n',hdr.qual.ac_max);
            
            fprintf(pfid,'\nAVG (summation) parmaters:\n');
            fprintf(pfid,'hdr.qual.mrdc_thr= %5.2f\n',hdr.qual.mrdc_thr);
            fprintf(pfid,'hdr.qual.mrac_thr= %5.2f\n',hdr.qual.mrac_thr);

            fprintf(pfid,'hdr.qual.dc_range= %5.2f\n',hdr.qual.dc_range);
            fprintf(pfid,'hdr.qual.ac_min= %5.2f\n',hdr.qual.ac_min);
            fprintf(pfid,'hdr.qual.delay_range= %5.2f\n',hdr.qual.delay_range);
            fprintf(pfid,'hdr.qual.mod_thr= %5.2f\n\n',hdr.qual.mod_thr);
            
            fprintf(pfid,'n_trials = %3.0f\n',n_trials);
            fprintf(pfid,'hdr.n_good.mrdc=     ');
            fprintf(pfid,'%4.0f',hdr.n_good.mrdc);
            fprintf(pfid,'\n');
            fprintf(pfid,'hdr.n_good.mrac=     ');
            fprintf(pfid,'%4.0f',hdr.n_good.mrac);
            fprintf(pfid,'\n');
            fprintf(pfid,'hdr.n_good.dc_range= ');
            fprintf(pfid,'%4.0f',hdr.n_good.dc_range);
            fprintf(pfid,'\n');
            fprintf(pfid,'hdr.n_good.ac=       ');
            fprintf(pfid,'%4.0f',hdr.n_good.ac);
            fprintf(pfid,'\n');
            fprintf(pfid,'hdr.n_good.ph=       ');
            fprintf(pfid,'%4.0f',hdr.n_good.ph);
            fprintf(pfid,'\n');
            fprintf(pfid,'hdr.n_good.mod=      ');
            fprintf(pfid,'%4.0f',hdr.n_good.mod);fprintf(pfid,'\n');
            fprintf(pfid,'hdr.n_good.tot=      ');
            fprintf(pfid,'%4.0f',hdr.n_good.tot);fprintf(pfid,'\n');
        end
        
        msg=strcat(msg,sprintf('; n-total =      '));
        msg=strcat(msg,sprintf('%5.0f',n2));
	
        % collect all n_good vals in one matrix and write to .xlsx file
        %
        fprintf('Collecting Params\n')
        params=[hdr.qual.ph_out_thr;hdr.qual.norm_mrdc_thr;hdr.qual.norm_mrac_thr;...
            hdr.qual.mrdc_thr;hdr.qual.mrac_thr;hdr.qual.dc_range;hdr.qual.ac_min;...
            hdr.qual.delay_range;hdr.qual.mod_thr];
        c_params={params};
        param_labels={'ph_out_thr';'norm_mrdc_thr';'norm_mrac_thr';...
            'avg_mrdc_thr';'avg_mrac_thr';'dc_range';'ac_min';...
            'delay_range';'mod_thr'};
        
        fprintf('Collecting n_goods\n')
        n_good=[1:hdr.n_chans;...
            good_norm_mrdc_flgs;good_norm_mrac_flgs;...
            hdr.n_good.mrdc;hdr.n_good.mrac;...
            hdr.n_good.dc_range;hdr.n_good.ac;...
            hdr.n_good.ph;hdr.n_good.mod;hdr.n_good.tot];
        c_good={n_good};
        n_good_labels={'Ch #';...
            'norm_mrdc';'norm_mrac';...
            'avg_mrdc';'avg_mrac';...
            'dc range';'ac min';...
            'delay range';'mod min';'total'};
        warning off MATLAB:xlswrite:AddSheet
        
        fprintf(['XLS Printing Params \n'])
        xlswrite([avg_PathName '\' xls_file_name],param_labels,['Block' num2str(i_block,'%02.0f')],'A1')
        xlswrite([avg_PathName '\' xls_file_name],params,['Block' num2str(i_block,'%02.0f')],'B1')
        
        fprintf(['XLS Printing Block# ' num2str(i_block) '\n'])
        xlswrite([avg_PathName '\' xls_file_name],n_good_labels,['Block' num2str(i_block,'%02.0f')],'A11')
        xlswrite([avg_PathName '\' xls_file_name],n_good,['Block' num2str(i_block,'%02.0f')],'B11:IW20')
        
	else % if various files exist or don't...
        if (in_file_flg < 0)
            msg=sprintf('Input file %s does not exist.\n',[in_path file_name]);
            fprintf(logid,'Input file %s does not exist.\n',[in_path file_name]);
        end
        if (evt_file_flg < 0)
            msg=sprintf('%s does not exist.\n',evt_file_name);
            fprintf(logid,'%s does not exist.\n',evt_file_name);
        end     
	end
	
else    % Out_file exists
    msg=sprintf('%s already exists!\n',[avg_PathName out_file_name]);
    fprintf(logid,'%s already exists!\n',[avg_PathName out_file_name]);
end
if print_flg==2;fclose(pfid); end
if logid ~= 1; fclose(logid); end